Projected northward shifts in eastern red‐backed salamanders due to changing climate

Abstract Many species' distributions are being impacted by the acceleration of climate change. Amphibians in particular serve numerous ecosystem functions and are useful indicators of environmental change. Understanding how their distributions have been impacted by climate change and will continue to be impacted is thus important to overall ecosystem health. Plethodon cinereus (Eastern Red‐Backed Salamander) is a widespread species of lungless salamander (Plethodontidae) that ranges across northeastern North America. To better understand future potential lungless salamander range shifts, we quantify environmental favorability, the likelihood of membership in a set of sites where environmental conditions are favorable for a species, for P. cinereus in multiple time periods, and examine shifts in the species' distribution. First, utilizing a large data set of georeferenced records, we assessed which bioclimatic variables were associated with environmental favorability in P. cinereus. We then used species distribution modeling for two time periods (1961–1980 and 2001–2020) to determine whether there was a regional shift in environmental favorability in the past 60 years. Models were then used to project future distributions under eight climate change scenarios to quantify potential range shifts. Shifts were assessed using fuzzy logic, avoiding thresholds that oversimplify model predictions into artificial binary outputs. We found that P. cinereus presence is strongly associated with environmental stability. There has been a substantial northward shift in environmental favorability for P. cinereus between 1961–1980 and 2001–2020. This shift is predicted to continue by 2070, with larger shifts under higher greenhouse gas emission scenarios. As climate change accelerates, it is differentially impacting species but has especially strong impacts on dispersal‐limited species. Our results show substantial northward shifts in climatic favorability in the last 60 years for P. cinereus, which are likely to be exacerbated by ongoing climate change. Since P. cinereus is dispersal‐limited, these models may imply local extirpations along the southern modern range with limited northward dispersal. Continued monitoring of amphibians in the field will reveal microclimatic effects associated with climate change and the accuracy of the model predictions presented here.


| INTRODUC TI ON
Anthropogenic climate change has been accelerating and will continue to accelerate over the coming decades, which will have profound impacts on both human and natural systems (Claeys et al., 2019;Nerem et al., 2018;Smith et al., 2015). Changes in a wide variety of variables such as temperature and precipitation (e.g., Stendel et al., 2021) have already impacted species ranges for many taxa across the world (e.g., Aragón et al., 2010;Báez et al., 2020;Hill et al., 2011). The capacity of a species to shift its range in response to environmental disturbance or climatic changes grants species the ability to persist over long time scales. For example, nonmigratory European butterflies were shown to shift northward between 35 and 240 km in the 20th century (Parmesan et al., 1999), and coastal pelagic fishes have shifted poleward by hundreds of kilometers in the past 20 years (Champion et al., 2021) in response to climate and environmental changes. However, given the acceleration of climatic shifts, as well as other factors such as land-use changes and urbanization, species may not be able to shift their ranges quickly enough and may suffer local extirpations. It is important to note that species are not all impacted equally (Estrada et al., 2015;Estrada, Morales-Castilla, et al., 2016) with dispersal-limited species being especially vulnerable (Enriquez-Urzelai, Bernardo, et al., 2019;Muñoz, Miller Hesed, et al., 2016). Therefore, a better understanding of current and future range shifts is critically important to conservation initiatives.
Among plethodontid salamanders, the eastern red-backed salamander (Plethodon cinereus) has been extensively studied behaviorally (Fleming et al., 2021;Hedrick et al., 2021;Hernández-Pacheco et al., 2019;Muñoz, Miller Hesed, et al., 2016;Sanchez et al., 2020) and has a large range that spans across the northeastern part of the United States and southeastern Canada (Petranka, 1998). Throughout their range, they are a prevalent, yet cryptic element in forest ecosystems (Welsh & Droege, 2001). This previous work has established the microclimatic preferences of P. cinereus, such as how air and soil temperature and precipitation impact detection rates of P. cinereus (Hernández-Pacheco et al., 2019), and how temperature impacts the species' growth rates and the onset of sexual maturity (Muñoz, Miller Hesed, et al., 2016). As a result of this wealth of research, P. cinereus has become a model species for understanding the impacts of microclimate on amphibian ecology, behavior, and dispersal (reviewed in Fisher-Reid et al., 2021). However, the effects of macroclimate on P. cinereus have been less fully explored, and the potential impacts of future climate on the species as a whole have yet to be studied via species distribution modeling, which may provide important insights into potential distributional changes over the next century (e.g., Zhang et al., 2020 in Chinese giant salamanders).
Dispersal-limited species are particularly vulnerable to climatic shifts, as rates of change are likely to be higher than the species' capacity to move to a more suitable habitat and expand its range (Estrada et al., 2015;Gibbons et al., 2000;Midgley et al., 2006;Muñoz, Miller Hesed, et al., 2016;Schneider & Root, 1998; but see Smith & Green, 2005). P. cinereus is dispersal limited and will likely need to respond to changing environmental conditions via alternative pathways such as phenotypic plasticity or vertical movement in the earth (rather than a horizontal movement to more suitable surface environments) (Muñoz, Miller Hesed, et al., 2016). Regardless, Plethodon species must come to the surface to feed and lay eggs (Petranka, 1998) and therefore the shifting surface conditions resulting from rapid climate change will likely have substantial impacts on P. cinereus' range and distribution.
Conserving species requires protecting the habitat within their distributions (Lomolino, 2004;Real et al., 2009). Species distribution models (SDMs) and ecological niche models are useful tools for understanding how a species' range may overlap with those of other species (Caravaggi et al., 2017;Real et al., 2009), how it has changed in the past (Gutiérrez-Rodríguez et al., 2017a, 2017bWarwick et al., 2021), and how it may change in the future under different climate projections (Báez et al., 2020;Pearson & Dawson, 2003;Real et al., 2010;Zhang et al., 2020), which may potentially lead to future species range expansions and contractions. SDMs have been used extensively to assess climate impacts on amphibian ranges (Barbosa & Real, 2012;Gutiérrez-Rodríguez et al., 2017a, 2017bReino et al., 2017;Sánchez-Montes et al., 2019), and it has been found that ectotherm distributions are more strongly and directly impacted by climactic changes, likely as a result of physiological constraints, compared with endotherms, which are more affected by indirect climate effects (Aragón et al., 2010).

K E Y W O R D S
conservation biogeography, ecological niche models, favorability, fuzzy logic, Plethodon, Plethodontidae, species distribution models

Biogeography
There are a wide variety of procedures to perform SDMs, one of them being the favorability function (Real et al., 2006), which assesses how the local presence probability differs from that expected by chance, due to the local environmental conditions, independently of the species' prevalence. Environmental favorability has the advantage of being commensurable and thus directly comparable across species and time periods (Acevedo & Real, 2012;Barbosa, 2015;Real et al., 2006), allowing the use of fuzzy logic to directly combine continuous predictions without applying thresholds. We aim to determine the macroclimatic factors that affect the range of P. cinereus, and how future climate change may impact P. cinereus by assessing past, present, and future environmental favorability across its range. These data will generate new insights into large-scale climatic impacts on P. cinereus and build on our current understanding of P. cinereus ecology. As P. cinereus is becoming a model system in amphibian ecology (Fisher-Reid et al., 2021), these data will be potentially helpful in better understanding range shift ecology in amphibians generally, especially lungless amphibians, in the face of rapidly accelerating climate change.
To address our goal of understanding how the distribution of P. cinereus may be impacted by future climate change, we first evaluated the bioclimatic factors related to the modern P. cinereus distribution and whether the climatic favorability for P. cinereus has shifted spatially between two time periods : 1961-1980 and 2001-2020. Then, utilizing two climate models (CCSM4, MIROC-ESM), each under four different greenhouse gas (GHG) emission scenarios, we used fuzzy logic to assess how changes in environmental favorability may impact the P. cinereus distribution by 2070. Specifically, we aim to (1) define the regions that are currently climatically favorable for P. cinereus across its range, based on SDMs built from data from the past 20 years (2001-2020) and determine which bioclimatic variables are influential in defining that range, (2) examine whether environmental favorability for P. cinereus has shifted in the last 60 years (with data from two time periods : 1961-1980 and 2001-2020), and (3) predict future range shifts and identify areas of projected expansions and contractions as a function of projected future climate. Based on previous data, we expect that P. cinereus' distribution will be strongly linked to precipitation-based climatic variables, as has been found more generally for amphibians (Aragón et al., 2010;Buckley & Jetz, 2007;. We further hypothesize that environmental favorability for P. cinereus will shift increasingly toward the north by 2070 under progressively higher GHG emission scenarios in both climate models. Finally, we hypothesize that the absolute amount of favorable area will decrease with increasingly higher GHG emissions.

| Data collection and climate variables
Plethodon cinereus occurrence data were downloaded from Global Biodiversity Information Facility (GBIF, 2021). This returned 106,424 spatially referenced records for P. cinereus as of September 2021 (https://doi.org/10.15468/ dl.d3j5ys). Data were then checked for errors so that only suitable occurrence data were used in analyses. This involved removing all duplicates, retaining only records of confirmed presences, and using the scrubr v. 0.4.0. package in R (Chamberlain, 2020), removing occurrence records where the latitude-longitude data were incomplete, not possible, or unlikely (e.g., at 0, 0 coordinates). Finally, given that records of P. cinereus should only occur in the United States and Canada, only records within those countries were accepted. Several records were found in western North America that were not removed earlier in the process and thus they were removed by not accepting any occurrence records west of a decimal longitude of −94 or south of a decimal latitude of 32. This resulted in 102,697 records. Data were not evenly distributed through time ( Figure S1a). The oldest record is from 1809, with the majority of records occurring before 1980 (79%) and after 2010 (only 7.1% of the records occurred between 1980 and 2010).
For evaluating shifts, we used subsets of the data for the periods 1961-1980-2020. The 1961-1980 data subset contained 58,621 records ( Figure S1b) and the 2001-2020 subset contained 12,521 records ( Figure S1c). Occurrence data were then formatted to be organized in a spatial data frame. A buffer of 350 km around the occurrence data was then established using the buffer function in the raster v. 3.5-29 R package (Hijmans, 2021) using the entire data set to establish environments that were accessible but where P. cinereus did not occur ( Figure S2). The buffered region was used as the study area. One method for delimiting a study area is to use geographical barriers that the species is not capable of crossing.
However, given the wide range of P. cinereus, we opted to have a buffered region that included inhospitable northern and southern climates as well as western climates where the species does not occur, but that are theoretically accessible to P. cinereus.
Environmental variables were downloaded from WorldClim (Hijmans et al., 2005) using the sdmpredictors v. 0.2.12 R package (Bosch, 2020) and then clipped to fit the buffered region outlined above. Each variable map had a spatial resolution of 0.08333 geo-  (Table 1). We excluded BIO3 (Isothermality), BIO14 (Precipitation of Driest Month), and BIO15 (Precipitation Seasonality) in our analyses due to a documented low correlation between present and future variables (Bedia et al., 2013). Occurrence data were then thinned with the gridRecords function of fuzzySim v. 4.3 package (Barbosa, 2015), to match the resolution of the climate variables.
Pixels that were not intersected by any presence points were used as unoccupied background.

| Bioclimatic modeling of P. cinereus between 1961-1980 and 2001-2020
We performed two approaches. We first ran simple models using all predictor variables and three different modeling techniques (glm, gam, and Maxent), to ensure that the choice of technique did not have a noticeable impact on the results. Then, we refined the models made with one of these techniques, namely glm, to conduct all the remaining analyses. First, generalized linear models (glm), generalized additive models (gam), and maximum entropy models (Maxent) were used to assess how bioclimatic variables were correlated with P. cinereus presence between 1961-1980 and 2001-2020. Two models were run for each modeling technique, one for each data subset. The predicted function of R was then used to generate presence probabilities across the study region. Correlations among the predictions of all three modeling techniques were high (Table S1, Figures S3 and   S4), thus it is unlikely that selecting a particular modeling method would strongly influence the predictions (Appendix S1). Therefore, additional analyses rely on glms, which produce presence probability and can be used to generate favorability (Acevedo & Real, 2012;Real et al., 2006).
To build a more refined glm, we used the multGLM function in the fuzzySim R package, using 20% of the data for model testing and 80% for model training for both data sets. This function removes highly correlated variables (among those pairs of variables with a Pearson correlation value above 0.8, it excludes the one with the least informative individual relationship with the distribution of the species). It also removes potentially irrelevant variables following a false discovery rate (Garcia, 2003), and then does a forward-backward stepwise selection of the remaining variables using an information criterion (AIC; Akaike, 1973;Burnham & Anderson, 2002). Nonsignificant variables after the stepwise procedure were removed using the model Trim function in fuzzySim (Crawley, 2007). This demonstrated which variables were significantly correlated with the P. cinereus spatial distribution data.
The output of the glm is a presence probability value. Probability depends both on the response of the species to the predictors and on the overall prevalence of the species (Cramer, 1999) (prevalence being the ratio between presences and the total number of cells). To remove the effect of prevalence from the model output, we applied the favorability function proposed by Real et al. (2006): where P is the probability value in a cell, n 1 is the total number of presences, and n 0 is the total number of absences in the data set.
The favorability function reflects the degree (between 0 and 1) to which the local probability values differ from that expected according to the species prevalence, where F = 0.5 corresponds to P = prevalence. Favorability values only reflect the response of the species to the predictors (Acevedo & Real, 2012) and (unlike probability) can be regarded as the degree of membership of the localities to the fuzzy set of sites with conditions that are favorable for the species (Acevedo & Real, 2012), which enables the easy application of fuzzy logic operations to distribution modeling (e.g., Robertson et al., 2004). Fuzzy logic operations expand the potential of the favorability function for comparison between models for different species, regions, and/or time periods (Estrada et al., 2008(Estrada et al., , 2011Real et al., 2010). Favorability was calculated using the Fav function in fuzzySim, using the proportion of presences included in the model (training data).
Following Reino et al. (2017), we assessed model performance in terms of discrimination capacity and calibration or reliability for both time periods. The threshMeasures function in the modEvA v.
3.5 package (Barbosa et al., 2013) was applied to assess the classification accuracy of the predicted distribution data based on the observed data. This was done using both prevalences as the classification threshold and then again using the maximum true skill statistic (TSS), which optimizes the sum of sensitivity and specificity.
Model evaluation metrics (correct classification rate, Sensitivity, Specificity, Precision, Cohen's kappa, and TSS) were calculated using each of these thresholds. Overall model discrimination capacity was assessed using the area under the curve (AUC) of the receiver operating characteristic (ROC) plot. The AUC gives the discrimination performance of the model over the entire range of prediction thresholds. Finally, the correlation between the observations and predictions was tested for significance, and Miller (Miller et al., 1991) calibration statistics were calculated (Wintle et al., 2005).

| Bioclimatic modeling of P. cinereus in the future
For future models, we utilized the WorldClim future layers in the sdmpredictors package (Hijmans et al., 2005)

| RE SULTS
There was substantial overlap in variables that were included in the glm models for the two time periods (Table 2). In the 1961-  (Swets, 1988 (Figures 1 and 2). There was also moderate expansion of favorability in the western part of the P. cinereus range.
Favorability contractions were most pronounced on the southeastern part of the P. cinereus range, especially in Virginia and the Carolinas. Fuzzy range change metrics suggested substantial gain relative to loss, so the expansion of the favorability range from the 1961-1980 to 2001-2020 time periods was larger than the associated contraction (Table S3).
Both the CCSM4 and MIROC-ESM climate models point to a potential northward shift in environmental favorability for P. cinereus  (Table S3). Fuzzy range change metrics for both climate models showed a similar pattern, whereby RCP 2.6 models had a relatively small range change while RCP 4.5 and 6.0 both had higher changes of a similar magnitude, and RCP 8.5 had substantial range gains and losses ( Figure 4).
Further, across all GHG pathways, the CCSM4 model suggested less change between the present and the future than the MIROC-ESM model.
Following previous studies, we focused on the RCP 8.5 models for further comparisons (Caravaggi et al., 2017).

| DISCUSS ION
Climate change is likely to cause major shifts in environmental suitability for a wide range of taxa, which will undoubtedly lead to range shifts and local extirpations (Pacifici et al., 2015;Van der Putten, 2012). Dispersal-limited species are likely to be more strongly impacted by these shifts because they cannot quickly respond to rapid environmental change. Given the wide geographical range of P. cinereus and the numerous investigations into its life history and behavior (reviewed in Fisher-Reid et al., 2021), as well as the importance of salamanders generally as environmental indicators (Buckley & Jetz, 2007;Feder & Burggren, 1992;Fleming et al., 2020;Welsh & Droege, 2001), this species is an excellent model system for utilizing species distribution models/ecological niche models to better understand the macroclimatic factors driving past, present, and future potential range shifts in dispersal-limited amphibians.
We found that a combination of precipitation and temperature variables define the regions that are currently climatically favorable for P. cinereus. In particular, regions with increased P. cinereus environmental favorability are more strongly associated with temperate, milder climates. Additionally, we found a northward shift in favorability between 1961-1980 and 2001-2020. Climate models indicate an additional northward shift by 2070. Given the low dispersal rates of P. cinereus and amphibians generally (Abellán & Svenning, 2014;Muñoz, Miller Hesed, et al., 2016;Recuero & García-París, 2011), this may suggest future extirpations.

F I G U R E 4
Fuzzy range change metrics for the two climate models (CCSM4, MIROC-ESM) under four GHG RCPs (2.6, 4.5, 6.0, and 8.5). Each plot shows gained favorability, lost favorability, stable favorability, stable absences, and balance (proportion for gained to lost favorability). Lost favorability generally increases with increasing greenhouse gas emissions.

| The impact of bioclimatic variables on P. cinereus' distribution
Due to their reliance on water for breeding, amphibians have been shown to be tied strongly to precipitation-based climate variables (Aragón et al., 2010;Buckley & Jetz, 2007;Duellman, 1999). Since most amphibians lay eggs in water, particularly in vernal pools, variation in precipitation can strongly impact amphibian reproductive success (Carey & Alexander, 2003). For example, variation in rainfall can change the number of eggs laid in a particular year (Caldwell, 1987). Although plethodontids are not tied to water for breeding, adult survivorship is also impacted by precipitation since they respire through their skin and have relatively high rates of water loss (Carey & Alexander, 2003;Shoemaker et al., 1992). Plethodontids lay their eggs in moist areas under rocks or fallen trees, with some species utilizing subterranean environments , and they spend time at the surface to feed (Petranka, 1998). As a result, P. cinereus must spend time at the surface and, given their physiological limitations, increased precipitation and soil moisture grant them the necessary conditions for surface activity (Gade & Peterman, 2019;Milanovich et al., 2010;Wilk et al., 2020). Further, experimental analyses have shown that increased temperatures reduce growth rates in P. cinereus (Muñoz, Miller Hesed, et al., 2016). Higher mean maximum July temperatures and lower precipitation in the driest month has been shown to be correlated with a 1.8% increase in P. cinereus body size between 1950-1970and 1980-2000(McCarthy et al., 2017. Therefore, changes to both precipitation and temperature variables are likely to strongly impact P. cinereus physiology. Our models suggest that favorable regions for P. cinereus in the present (2001-2020) are centered on climatic stability (Table 2). P. cinereus presence was correlated with higher precipitation in the driest quarter (BIO17), lower precipitation in the wettest quarter (BIO16), lower annual temperature range (BIO7), lower temperature in the driest quarter (BIO9), lower precipitation in the warmest quarter (BIO18), and higher mean temperature in the wettest quarter (BIO8).
These factors combined do not suggest that higher precipitation implies the higher chance of presence, but rather that moderate, mild climates lead to higher presence probability and more climatically favorable regions for P. cinereus. This may relate to a wider pattern of niche conservatism in climatic tolerances found in plethodontids in Appalachia, whereby the highest species diversity is found in intermediate-elevation habitats (Kozak & Wiens, 2010). Kozak and Wiens (2010) suggest that plethodontids are specialized to a certain set of environmental conditions leading to a build-up of species in regions mimicking their ancestral environments due to niche conservatism. Climate change is leading to a greater frequency and severity of both precipitation and temperature extremes (Stott, 2016;Zhang et al., 2013). Therefore, this has strong implications for P. cinereus distributions in the future. Since P. cinereus spends most of its life underground, soil temperature 30 cm below the surface (rather than air temperature or precipitation) serves as the best metric for predicting detection rates (Sanchez et al., 2020). P. cinereus are likely able to modulate their position in the soil column so that they are at optimum environmental conditions (Muñoz, Miller Hesed, et al., 2016). However, P. cinereus does need to come to the surface to feed and mate (Houck & Verrell, 1993;Tilley & Bernardo, 1993). More variable macroclimate conditions will likely have strong effects on favorable surface conditions for salamanders, and in turn impact their behavior and physiology (Muñoz, Miller Hesed, et al., 2016). Warmer temperatures have been shown to negatively impact autumn growth in P. cinereus, increasing the time until individuals are sexually mature (Muñoz, Miller Hesed, et al., 2016). Surface activity is impacted partly by macroclimatic variables, but also by microclimatic fluctuations. For example, canopies buffer against microclimatic variability (Chen et al., 1999), and therefore, canopy loss may have strong local effects on salamanders (Demaynadier & Hunter, 1998;Harpole & Haas, 1999).

|
Many salamander species, including P. cinereus, have generally small home ranges and low dispersal capabilities (Abellán & Svenning, 2014;Lunghi & Bruni, 2018;Muñoz, Miller Hesed, et al., 2016;Recuero & García-París, 2011). Indeed, previous work has shown that P. cinereus individuals are usually found within a meter of their original capture site (Gergits & Jaeger, 1990 (Kearney & Porter, 2009;Thomas et al., 2004) and low dispersal capabilities may lead to results that conflict with distribution model predictions (Pearson & Dawson, 2003). As a result, northward shifts in environmental favorability for P. cinereus do not necessarily imply northward shifts in the species' distribution range. Our data may suggest extirpation (if no local adaptation occurs) in the southern regions as they become less environmentally favorable, without the coincident northward range expansion implied by our models (Figure 3). Similarly, these models are correlative and do not include the fact that P. cinereus may modify its behavior or be physiologically plastic to persist in environmentally unfavorable conditions at the southern edge of its range (Newman et al., 2022;Riddell & Sears, 2020). As a result of these factors, the shifts in favorability should be best viewed as potentially suitable areas and future work may improve our projections with additional analyses.
On the contrary, other amphibian distribution modeling work has found that Lissotriton helveticus was able to shift northward during rapid warming in the Holocene despite low dispersal capabilities (Recuero & García-París, 2011). Other amphibians such as Hyla molleri were not strongly impacted by climate fluctuations since the Last Interglacial, since it has a relatively high cold tolerance (Sánchez-Montes et al., 2019). Further, recent work has shown that the majority of P. cinereus genetic diversity is in the southern part of its modern range (Radomski et al., 2020). Three quarters of the modern P. cinereus range was covered by glaciers during the last glacial maximum (LGM) and as such, rapid range expansion after the LGM is postulated (Radomski et al., 2020). Plethodon albagula was shown to have relatively low gene flow through favorable habitats, suggesting that plethodontids may move quickly and directly through inhospitable habitats while they move in a slower, more exploratory fashion in favorable habitats (Peterman et al., 2014), potentially increasing dispersal capabilities when pressed. Thus, P. cinereus and other amphibians may have higher dispersal capabilities than currently realized. Further, the large population size of P. cinereus may enhance the likelihood of persistence and future recolonization from climate refugia following rapid climatic shifts (Recuero & García-París, 2011). It is noteworthy that climate is often found to be less impactful on species than direct anthropogenic changes (Gutiérrez-Rodríguez et al., 2017b;Hampe & Jump, 2011). Matching SDM and phylogeographic data with dispersal simulations (e.g., Pearson & Dawson, 2003;Peterson et al., 2001) would be a fruitful avenue for better understanding how dispersal-limited species will react to future shifts in climate favorability.

ACK N OWLED G M ENTS
We thank the members of the Salamander Population and Nacionales of Spain through project 2745/2021. We also thank three anonymous reviewers and the editor for greatly improving the quality of the manuscript.

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declare that they have no conflict of interest.